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Addendum  to  Final  Report 


AFOSR  Grant  FA9550-07-1-0047 
Human  Supervision  of  Time  Critical  Control  Systems 


Jordan  Cannon 
Pavlo  Krokhmal 
Robert  Murphey 
Panos  Pardalos 


Abstract.  This  addendum  reports  findings  of  the  Pis  in  the  course  of  the  project  “Human 
Superv  ision  of  Time-Critical  Control  Systems”,  which,  at  the  moment  of  the  report,  have 
not  been  published  or  submitted  for  publication.  In  particular,  it  presents  a  robust 
technique  for  detecting  temporal  changes  in  multidimensional  time  scries  data 
represented  by  phychophysiological  measurements.  The  proposed  approach  relics  on 
detecting  trends  in  an  appropriate  statistics  obtained  using  independent  component 
analysis  of  the  data.  In  comparison  to  the  previous  results  obtained  by  the  Pis  in  the  scope 
of  this  project,  the  new  method  performs  much  more  robustly  in  terms  of  between-subject 
and  between-trial  differences.  The  proposed  algorithm  is  amenable  to  efficient 
implementation  compliant  with  on-line  performance  requirements. 
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1.  Introduction 

The  increasing  complexity  and  sophistication  of  computing,  sensing,  and 
communication  technologies  paves  the  way  for  proliferation  of  unmanned  autonomous 
systems  and  platforms,  which  will  replace  and/or  assist  humans  in  hazardous  or  resource 
consuming  missions.  Yet,  the  common  consensus  is  that  despite  the  ever  growing  need 
for  increased  autonomy  of  various  unmanned  systems  and  vehicles,  the  human 
supervision  of  sueh  systems  is  indispensable  and  eritieal  for  mission  suecess. 

On  the  other  hand,  the  large  amounts  of  information  produced  by  these  eomplex 
systems  can  place  a  high  demand  on  a  human  operator’s  cognitive  load,  potentially 
overwhelming  him/her  and  leading  to  degraded  performance.  Thus,  the  ability  of  an 
automated  control  system  to  estimate  the  current  functional  state  of  a  human 
operator/supervisor  and  supply  information  that  is  conducive  to  human  decision-making 
at  the  given  cognitive  level  is  essential  to  robust  and  sueeessful  system  performance. 

The  objective  of  the  present  endeavor  is  to  advance  methods  for  real-time 
detection  of  changes  in  human  operator’s  functional  state  that  meet  real-time  or  online 
requirements,  and  which  can  be  utilized  in  elosed-loop  autonomous  control  systems  with 
human  supervision.  Generally,  operator  functional  state  (OFS)  can  be  defined  as  the 
momentary  ability  of  an  operator  to  meet  task  demands  with  their  cognitive  and 
physiological  resources.  In  the  context  of  this  work,  the  OFS  is  associated  with  cognitive 
load  experienced  by  the  operator;  with  this  eaveat  in  mind,  we  will  use  both  terms 
interchangeably. 

With  the  ability  to  measure  and  deteet  changes  in  OFS  in  real-time,  a  elosed-loop 
system  between  the  operator  and  machine  eould  optimize  OFS  through  the  dynamie 
allocation  of  tasks.  For  instance,  if  the  system  detects  the  operator  is  in  cognitive 
overload,  it  ean  automate  eertain  tasks  allowing  them  to  better  foeus  on  salient 
information.  Conversely,  if  the  system  detects  under-vigilance,  it  ean  allocate  tasks  back 
to  the  manual  control  of  the  operator.  In  essenee,  this  system  operates  to  “dynamically 
match  task  demands  to  [an]  operator’s  momentary  cognitive  state”,  thereby  achieving 
optimal  OFS  (Wilson,  Russell,  2007). 

OFS  is  eommonly  measured  indirectly,  like  using  overt  performance  metrics  on 
tasks;  if  performance  is  declining,  a  low  OFS  is  assumed.  Another  indirect  measure  is  the 
subjective  estimate  of  mental  workload,  where  an  operator  narrates  his/her  perceived 
functional  state  while  performing  tasks  (Wilson,  Russell,  2007).  Unfortunately,  indirect 
measures  of  OFS  are  often  infeasible  in  operational  settings;  performance  metries  are 
difficult  to  construct  for  highly-automated  eomplex  systems,  and  subjective  workload 
estimates  are  often  inaccurate  and  intrusive  (Wilson,  Russell,  2007;  Prinzel  et  al.,  2000; 
Smith  et  al.,  2001). 

OFS  can  be  more  directly  measured  via  psychophysiological  signals  sueh  as 
eleetroeneephalogram  (EEG)  and  eleetrooculography  (EOG).  Current  research  has 
demonstrated  these  signals’  ability  to  respond  to  changing  cognitive  load  and  to  measure 
OFS  (Wilson,  Fisher,  1991;  Wilson,  Fisher,  1995;  Gevins  et  al.,  1997;  Gevins  et  al., 
1998;  Byrne,  Parasuraman,  1996).  Moreover,  psyehophysiologieal  signals  arc 
continuously  available  and  ean  be  obtained  in  a  non-intrusive  manner,  pre-requisite  for 
their  use  in  operational  environments. 

The  objective  of  this  study  is  to  advanee  schemes  whieh  deteet  ehanges  in  OFS  by 
monitoring  psychophysiological  signals  in  real-time.  Reviews  on  similar  methods  ean  be 
found  in,  e.g.,  Wilson  and  Russell  (2003a)  and  Wilson  and  Russell  (2007).  Many  of  these 
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methods  employ  pattern  recognition  to  classify  mental  workload  into  one  of  several 
discrete  categories.  For  instance,  given  an  experiment  with  easy,  medium  and  hard  tasks, 
and  assuming  the  tasks  induee  varying  degrees  of  mental  workload  on  a  subjeet,  these 
methods  classify  whieh  task  is  being  performed  for  each  epoeh  of  psyehophysiologieal 
data.  The  most  eommon  elassifiers  are  artificial  neural  networks  (ANN)  and  multivariate 
statistical  techniques  sueh  as  stepwise  discriminant  analysis  (SWDA).  ANNs  have  proved 
especially  effective  at  classifying  OFS  as  they  aeeount  for  the  non-linear  and  higher  order 
relationships  often  present  in  EEG/EOG  data;  they  routinely  aehieve  classification 
aeeuraey  greater  than  80%. 

However,  the  discrete  output  of  these  classification  schemes  is  not  eondueive  to 
real-time  ehange  detection.  They  accurately  classify  OFS,  but  they  do  not  indicate  when 
OFS  has  changed;  the  ehange  points  remain  ambiguous  and  left  to  subjective 
interpretation.  Thus,  the  present  study  introduces  several  online  algorithms  whieh 
objectively  determine  ehange  in  OFS  via  real-time  psyehophysiologieal  signals. 

The  following  sections  describe  the  dataset  evaluated,  discuss  the  statistical 
properties  of  psyehophysiologieal  signals,  and  detail  the  various  algorithms  whieh  utilize 
these  signals  to  deteet  real-time  ehange  in  OFS.  For  each  algorithm,  results  are  presented 
regarding  their  effieaey  and  a  diseussion  is  provided.  Finally,  the  study  is  eoneluded  with 
a  review  and  comparison  of  each  method. 

2.  Data  collection,  processing,  and  analysis 

The  dataset  utilized  in  the  following  analyses  originated  from  experiments 
eondueted  at  Wright-Patterson  AFB  in  2008.  Data  was  available  for  three  subjects.  Eaeh 
subjeet  performed  two  fourteen-minute  trials  with  unmanned  aerial  vehicle  (UAVj  tasks 
presented  at  three  levels:  low,  medium,  and  high,  denoted  as  LL,  ML,  and  HL, 
respectively.  The  LL  was  the  baseline  state  and  subsequently  encompassed  most  of  eaeh 
trial.  The  ML  and  HL  were  presented  four  times  eaeh,  in  a  balanced  order,  lasting 
approximately  20  seconds  eaeh  time.  Eaeh  trial  began  in  the  LL  and  after  every  ML  or 
HL  was  presented,  it  returned  to  the  LL.  The  experimental  design  assumed  that  varying 
task  load  indueed  corresponding  levels  of  cognitive  load  on  the  subjeet.  Throughout  this 
report,  cognitive  load  is  assumed  to  be  a  proxy  measure  of  OFS;  thus  detecting  a  ehange 
in  task  load,  is  synonymous  with  detecting  a  ehange  in  OFS. 

The  tasks  involved  monitoring  four  UAVs  executing  a  bombing  mission.  During 
the  ML  and  HL,  subjects  performed  a  visual  search  of  a  radar  image  uploaded  at 
designated  waypoints  by  the  UAVs.  The  subjects  panned  the  radar  image,  loeated 
potential  targets,  and  marked  six  of  them  for  bombing,  according  to  pre-determined 
priorities.  This  had  to  be  accomplished  within  the  20  second  timeframe.  The  HL  was 
more  difficult  than  the  ML  as  its  radar  image  contained  more  distraeter  targets,  in 
addition  to  imposing  more  complex  targeting  priorities.  The  baseline  LL  condition  simply 
required  the  subjeet  to  monitor  the  UAVs  flight  paths  until  eaeh  waypoint,  i.e.  ML  or  HL, 
was  encountered.  Overall,  the  tasks  were  very  visual  in  nature  and  were  expected  to 
engage  the  visual  processing  eenters  of  the  brain. 

The  data  eolleeted  for  eaeh  trial  consisted  of  eight  physiological  channels  of  EEG, 
EOG,  and  eleetroeardiogram  (ECG)  recorded  at  a  sampling  frequency  of  200  Hz.  The 
EEG  channels  were  recorded  from  five  electrodes:  F7,  Fz,  Pz,  T5,  and  O2.  These  electrodes 
were  affixed  to  the  subject’s  sealp  aeeording  to  the  10/20  International  electrode  system 
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shown  in  Figure  2.1.  Vertical  and  horizontal  EOG  data,  termed  VEOG  and  HEOG 
respectively,  were  collected  for  two  purposes:  primarily,  as  a  measure  of  cognitive  load, 
and  secondly,  to  eliminate  blink  artifacts  in  the  EEG  signals.  Finally,  one  channel  of  ECG 
was  collected  to  measure  heart  rate.  Only  the  EEG  and  EOG  data  were  used  in  the 
following  analyses. 

From  this  point  forward,  each  trial  will  be  denoted  by  a  letter  followed  by  a 
number;  the  letter  represents  the  subject  A,  E,  or  F,  and  the  number  identifies  whether  it 
was  the  subject’s  first  or  second  trial.  For  instance,  A01  denotes  data  from  the  first  trial  of 
subject  A. 


Front 


Back 


Figure  2.1  EEG  electrode  diagram 


2.1.  Preprocessing  the  data 

Blink  artifacts  contaminate  EEG  signals  when  the  electrical  activity  from  a 
subject’s  eye  is  recorded  by  the  EEG  electrodes  affixed  to  their  scalp.  This  interference 
greatly  distorts  the  EEG  data  and  negatively  affects  subsequent  analyses.  Many  filters 
have  been  developed  to  remove  blink  contamination,  of  interest  here  however,  are  filters 
that  remove  the  artifacts  online,  rather  than  in  a  post-processing  scheme.  One  form  of 
online  artifact  removal  is  adaptive  filtering.  The  adaptive  filter  incorporates  the  VEOG 
and  HEOG  signals  as  reference  inputs  to  de-contaminate  an  EEG  signal,  s(i),  for  every 
time  moment  i.  The  artifact-free  signal,  e(i),  results  by 


e(0  =  s(i)  -  rv(i)  -  rh(i) 


(2.1) 


where, 
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M 


fv(i)  =  ^  hv(m)rv(n  +  1  -m) 


m= 1 


M 


*h(0  =  ^  hh(m)rh(n  +  1  -  m) 


m=l 


(2.2) 


are  filtered  VEOG  and  HEOG  reference  signals,  respectively.  The  hv{m )  and  hh(m) 
terms  represent  finite  impulse  response  (FIR)  filters  of  length  M,  which  are  updated  for 
every  time  period  i,  to  filter  the  raw  VEOG  and  HEOG  signals. 

Updating  the  FIRs  is  accomplished  through  a  recursive  least-squares  (RLS) 
algorithm  presented  in  Table  2.1.  The  underbars  denote  column  vectors  and  /?(i)  is  a 
matrix.  In  the  present  analysis,  the  forgetting  factor,  X,  was  set  to  .9999  and  M  was  set  to 
a  length  of  three.  See  He,  Wilson,  and  Russell  (2007)  for  complete  details. 

Figure  2.2  (a)  displays  an  EEG  signal  recorded  at  the  F7  electrode,  and  (b) 
displays  the  VEOG  signal  recorded  during  the  same  period  of  time.  Blink  contamination 
is  clearly  present  in  (a),  where  sharp  peaks  indicate  eye  blinks.  Notice  these  same  blinks 
are  reflected  in  (b),  the  VEOG  signal,  whose  intent  is  to  detect  eye  activity.  The  adaptive 
filter  uses  the  VEOG  signal  to  identify  the  blinks  and  to  remove  them  from  affected  EEG 
signals,  like  F7,  where  blink  activity  is  considered  interference.  Figure  2.3  depicts  the 
contaminated  F7  signal  in  (a),  but  this  time  with  the  corresponding  filtered  F7  signal  in 
(b),  after  adaptive  filtering.  Notice  that  the  sharp  peaks  characteristic  of  eye  blinks  arc  no 
longer  present. 

Once  eye  blink  artifacts  were  removed,  the  psychophysiological  signals  were 
subjected  to  Discrete  Fourier  Transform  (DFT)  for  every  epoch,  usually  three  to  five 
seconds.  The  DFT  transforms  signals  from  the  time  domain  to  the  frequency  domain.  To 
accomplish  this,  the  DFT  assumes  that  the  time  domain  signal  is  a  sum  of  many 
sinusoids;  this  assumption  is  generally  deemed  appropriate  for  psychophysiological 
signals.  The  DFT  and  its  inverse  are  computed  by 


JV(fl)  = 


00 

^  x[n\e~>na 

n=-oo 


i  r2n 

x[n]  =  —  I  X(CL)e’na  dCl 

2nJ0 


(2.3) 


where  fl  is  the  discrete-frequency  variable  and  x[n]  is  a  discrete  time  series,  in  this  case 
obtained  by  sampling  a  continuous  EEG^EOG  signal  (Phillips  et  al.,  2007). 

Utilizing  the  DFT,  a  frequency  spectrum  was  created  for  each  epoch  where  the 
powers  of  particular  wavebands,  e.g.  8-12  Hz,  known  as  alpha,  beta,  etc.,  were 
computed.  Waveband  power  was  computed  as  the  average  of  each  frequency’s  power,  for 
frequencies  falling  within  the  waveband.  The  powers  of  particular  EEGYEOG  wavebands 
have  been  shown  to  correlate  with  changes  in  cognitive  load  (see  Table  2.2).  For 
instance,  alpha  waveband  power,  8-12  Hz,  increases  with  relaxation  while  theta 
waveband  power,  5-8  Hz,  generally  decreases  with  relaxation  (Smith  et  ah,  2001). 
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Table  2.1  Adaptive  filtering  algorithm 


1.  Initialize: 
1.1.  H( 


H_(i  —  1)  =  0  where  H_  =  j*  and  and  arc  vectors  of  the  filter 


coefficients,  hv(m )  and  hh(rn),  respectively 

1.2.  [R(i  —  l)]-1  =  //< 7  where  /  is  a  2Mx  2 M  identity  matrix  and  a  =  0.01 

1.3.  M  and  \  are  user-defined 

1.4.  i  =  M 

1.5.  n  is  the  length  of  the  signal  to  be  filtered 
2.  Calculate  /£(i): 


[R(i  -  DJ-VCO 


A  +  r(i)T[R(i  -  l)]_1r(i) 


T*  ( 

where  r(i)  =  —  where  rv(i)  and  rh(i)  are  vectors  of  the  VEOG  and  HEOG  signal, 

—  — 

respectively: 


rv(i)  =  [rv(i),rv{i  -  1), ...  ,rv(i  +  1  -  M)]r 


rh(0  =  [rn(0,rh(i  -  1),  ...,rh(i  +  1  -  M)]T 


3.  Calculate  e 


4.  Calculate //(i): 


5.  Update  [/?(0] 


[Rii)]-1  =  A-1[/?(t  —  1)]_1  -A-1ff(i)r(Or[«(i-  l)]"1 


6.  Calculate  e(i)  : 


e(Q  =  s(Q-r(QTH(Q 
Note:  This  is  the  vector  equivalent  of  (2. 1) 


7.  i  =  /+l 


8.  While  i  <  n  repeat  steps  2  through  8 


100 


7 
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Figure  2.3  (a)  Un-filtered  and  (b)  filtered  F7  signals 
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Tabic  2.2.  EEG  wavebands 


Waveband 

Frequency  (Hz) 

Interpretation 

Delta 

1-4 

Slow  wave  sleep 

Theta 

4-8 

Arousal 

Alpha 

8-12 

Relaxation 

Beta 

12-30 

Active,  alert,  working 

Gamma 

30+ 

Cognitive  and/or  motor  function 

From  this  point  forward,  the  waveband  powers  computed  from 
psychophysiological  signals  recorded  at  particular  electrodes  arc  referred  to  as  features.  It 
is  possible  to  compute  hundreds  of  features  for  each  epoch,  as  there  are  many  different 
combinations  of  wavebands  and  electrodes.  Methods  which  measure  and  detect  real-time 
changes  in  OFS  must  determine  which  of  these  features  are  indicators  of  cognitive  load, 
and  which  are  noise. 


2.3.  Statistical  Analysis  of  the  Post-Processed  Data 

Empirical  analysis  was  conducted  on  several  psychophysiological  features  to 
characterize  their  properties  and  behavior  with  respect  to  subject  and  task  load.  For  this 
analysis,  the  features  were  computed  in  three-second  epochs.  It  is  common  to  find  one- 
second  epochs  in  the  literature,  however,  they  often  exhibit  highly  variable  and  erratic 
behavior.  Moreover,  short  epochs  induce  autocorrelation  in  the  features.  In  contrast, 
longer  epochs  are  more  stable,  but  if  they  are  too  long,  they  will  be  infeasible  in 
operational  settings  where  change  occurs  instantaneously.  Three-second  epochs  were 
found  to  achieve  the  best  balance  between  short  and  long  epochs. 

The  features  chosen  for  this  analysis  are  traditional  features  used  in  OFS 
classification  methods.  The  wavebands  of  these  features  arc:  delta,  2-4  Hz,  theta,  5-8  Hz, 
alpha,  9-13  Hz,  beta,  14-32  Hz,  and  gamma,  33-43  Hz.  Features  from  these  wavebands 
have  shown  to  correlate  strongly  with  cognitive  load,  thus  they  are  often  used  to  classify 
OFS. 

Figure  2.4  displays  the  distributions  of  theta,  alpha,  and  beta  waveband  powers 
computed  across  all  electrodes  for  every  epoch.  The  data  from  every  trial  were  combined 
to  form  the  distributions,  after  the  data  were  standardized  to  zero  mean  and  unit  variance 
with  respect  to  each  trial.  As  shown,  the  distributions  are  non-normal,  displaying  a  strong 
positive  skew.  The  beta  waveband  power  is  especially  skewed  and  nearly  resembles  a 
Poisson  distribution.  These  distributions  were  subjected  to  statistical  tests  for  non¬ 
normality  and  each  was  strongly  significant  at  the  .05  level,  indicating  that  features  from 
these  wavebands  originate  from  non-normal  distributions. 

Recall  that  all  seven  electrodes  record  the  same  five  wavebands  for  every  epoch. 
The  reason  for  this  seemingly  “redundant”  data  is  that  brain  waves  behave  differently 
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when  emitted  from  different  regions  of  the  brain.  For  instance,  alpha  power  is  greatest  in 
parietal  regions,  e.g.  Pz,  whereas  theta  power  dominates  the  frontal  region,  e.g.  Fz.  While 
differences  exist  between  topographical  locations,  there  still  remains  significant 
correlation,  or  dependency,  between  identical  wavebands  collected  at  different 
electrodes;  hence  not  every  feature  records  unique  information. 

To  illustrate  the  correlation  between  features,  Table  2.3  contains  the  correlation 
coefficients  of  features  recorded  at  VEOG  of  F01.  Notice  that  every  feature  was 
positively  correlated  with  an  average  correlation  coefficient  of  .819.  In  this  sense,  one 
feature  from  FOl’s  VEOG  electrode  could  virtually  represent  all  the  information  recorded 
there.  The  understanding  of  this  phenomenon  implores  selecting  a  variety  of  features, 
from  different  electrodes,  so  they  each  contain  unique  information  about  the  subject’s 
cognitive  state. 


Figure  2.4  Histograms  of  theta,  alpha,  and  beta  wavebands  collapsed  across  all  subjects 


Many  of  the  features  exhibited  another  form  of  correlation,  called  autocorrelation. 
Autocorrelation  occurs  when  an  observation,  ,  recorded  at  time  i+k,  is  dependent 
upon  a  previous  observation,  ,  where  k  is  the  time  lag.  Thus,  if  an  autocorrclated  feature 
is  observed  at  time  i,  it  will  contain  information  about  its  value  k  time  moments  in  the 
future  (Montgomery  et  al.,  2008).  Autocorrelation  is  identified  by  evaluating  the  sample 
autocorrelation  function  for  a  finite  time  series,  ,  defined  by 


(2.4) 
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where  ck  and  c0  are  computed  from 

1  n~k 

ck=-  V  (x,-  -x)(xi+k  -x)  fork  =  0,1,2,  ...,K  (2.5) 

n  L—t 

i=l 

Beta  and  gamma  features  consistently  exhibited  autocorrelation.  Figure  2.5 
displays  the  sample  autocorrelation  function  for  a  beta  feature  measured  at  T5  of  E01.  At 
lag  1,  the  autocorrelation  is  the  strongest  and  then  exponentially  decreases  before 
oscillating  around  0.  This  behavior  is  characteristic  of  a  first-order,  autoregressive  time 
series  (AR[1]). 

A  survey  of  feature  autocorrelation  was  conducted  across  subjects.  The  results  arc 
summarized  in  Table  2.4.  For  each  waveband,  a  percentage  is  displayed  representing  the 
proportion  of  electrodes  where  the  feature  exhibited  significant  autocorrelation.  For 
instance,  the  alpha  waveband  for  A01  generated  features  which  were  autocorrelated  in 
three  of  the  seven  electrodes,  i.e.,  42.86%.  As  shown,  the  wavebands  which  generated  the 
least  autocorrelated  features  across  subjects  were  theta  and  alpha.  In  contrast,  beta  and 
gamma  wavebands  gave  rise  to  features  that  were  autocorrelated  regardless  of  subject  or 
electrode. 


Table  2.3  Inter-node  correlation  measured  at  VEOG  of  F01 


Delta 

Theta 

Alpha 

Beta 

Gamma 

Delta 

1.000 

0.922 

0.939 

0.852 

0.582 

Theta 

0.922 

1.000 

0.949 

0.864 

0.655 

Alpha 

0.939 

0.949 

1.000 

0.914 

0.663 

Beta 

0.852 

0.864 

0.914 

1.000 

0.853 

Gamma 

0.582 

0.655 

0.663 

0.853 

1.000 

It  is  also  important  to  characterize  how  the  properties  of  psychophysiological 
features  change  with  time.  More  specifically,  it  is  necessary  to  investigate  if  these 
features  exhibit  stationary  behavior.  Stationarity  implies  a  statistical  equilibrium,  where 
the  properties  of  a  time  series,  such  as  its  autocorrelation  or  probability  distribution,  arc 
stable  over  time.  A  time  series  is  strictly  stationary  when  the  joint  probability  distribution 
of  observations  xitxi+1,  ...,xi+n  is  the  same  as  the  joint  distribution  of 
xi+klxi+k+1, ... ,x j+fc+n-  However,  it  is  often  sufficient  to  classify  a  time  series  as 
stationary  as  long  as  it  varies  around  a  fixed  mean  (Montgomery  et  al.,  2008). 
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Figure  2.5  Sample  autocorrelation  function  for  T5  beta  of  E01 
Note:  the  red  lines  denote  .05  significance  levels 


Table  2.4  Results  of  autocorrelation  survey 


A01  E01  F01 

Average 

Delta 

0.00%  57.14%  71.43% 

14.29%  28.57%  42.86% 

42.86%  28.57%  42.86% 

100.00%  100.00%  85.71% 

100.00%  100.00%  100.00% 

42.86% 

28.57% 

38.10% 

95.24% 

100.00% 

Theta 

Alpha 

Beta 

Gamma 

By  the  latter  definition,  non-stationary  features  will  exhibit  some  type  of  trend,  for 
instance,  their  mean  might  increase  with  time.  These  trends  can  be  determined  by  plotting 
a  feature’s  autocorrelation  function.  If  the  plot  displays  slowly-decreasing  autocorrelation 
with  increasing  lag,  a  feature  is  classified  as  non-stationary.  When  this  analysis  was 
conducted  on  a  sample  of  features  from  each  waveband,  none  of  the  features  exhibited 
the  symptoms  of  non-stationarity.  Instead,  they  each  varied  around  a  fixed  mean  as  in  the 
example  shown  in  Figure  2.6. 

Finally,  an  empirical  analysis  was  conducted  to  characterize  the  relationship 
between  features  and  varying  task  load;  these  results  are  in  Appendix  A.  Although  no 
casual  relationship  between  psychophysiological  features  and  task  load  are  firmly 
established,  many  studies  have  asserted  strong  associations,  such  as  increasing  alpha 
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power  with  decreasing  task  load.  The  features  that  arc  most  valuable  to  OFS  classifiers, 
are  those  that  monotonically  increase  or  decrease  with  changing  task  load.  If  this 
relationship  is  present,  then  these  features  can  provide  information  on  a  subject’s 
cognitive  load.  Moreover,  it  would  be  ideal  if  particular  features  had  the  same  correlation 
to  task  load  regardless  of  the  subject;  however,  the  results  of  this  empirical  analysis  found 
no  such  “universal”  feature.  Instead,  the  features  of  each  subject  exhibited  unique 
behavior  with  respect  to  task  load. 

Despite  the  inability  of  features  to  generalize  across  subjects,  there  existed  at  least 
one  feature  for  each  subject  which  monotonically  increased  or  decreased  with  task  load. 
Thus  by  selecting  features  tailored  to  each  subject,  it  is  possible  to  gain  reliable 
information  on  OFS.  In  fact,  nearly  all  OFS  classifiers  adhere  to  this  “subject-specific” 
approach  (Wilson,  Russell,  2003a;  Smith  et  al.,  2001). 

2r 


o 


J. 


± 


SnrfeNuvter 


Figure  2.6  Stationary  plot  of  Fz  theta  for  A01 


3.  Independent  Component  Analysis  of  Multidimensional 
Psychophysiological  Time  Series 

Independent  component  analysis  (ICA)  is  a  method  used  to  extract  the  underlying 
components  of  signals.  ICA  assumes  that  physical  processes,  like  the  brain,  arc 
comprised  of  distinct  operators  which  emit  signals  independent  of  each  other.  When  these 
signals  are  recorded  by  sensors,  they  become  “mixed”  and  indistinguishable.  ICA  is  a 
method  for  separating  these  mixed  signals  into  their  underlying  components  (Stone, 
2002). 
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For  instance,  two  EEG  electrodes  record  brain  activity  at  different  locations  on 
the  scalp.  The  signals  recorded  from  these  electrodes,  denoted  as  *i(0  and  x2(i),  are 
assumed  to  be  a  mixture  of  components,  s1(i)  and  s2(i),  originating  from  independent 
operators  in  the  brain.  xt(i)  and  x2(i )  are  expressed  as  linear  combinations  of  the  source 
components, 


*i(0  =  %iSi(i)  +  a12s2(0 

*2(0  =  «2iSt  (i)  +  a22s2(i) 


(3.6) 


The  objective  of  ICA  is  to  reveal  the  source  components  using  only  the  recorded 
EEG  signals,  xx (i)  and  x2(i).  As  shown  in  (3.6),  this  requires  determining  the 
coefficients  a^.  This  is  a  difficult  task,  however,  it  is  made  feasible  by  assuming  the 
source  components  are  independent  of  each  other,  hence  “independent  components” 
(Hyvarinen,  Oja,  2000). 

ICA  is  more  generally  defined  by  considering  the  system  of  equations  in  (3.6)  in 
vector-matrix  notation, 


x(i )  =  As(i )  (3-7) 

where  A  is  the  mixing  matrix,  the  recorded  psychophysiological  signals  are  in  vector 
x(i),  and  the  independent  components  arc  in  vector  s(i),  for  each  epoch  i.  Once  A  is 
determined,  its  inverse,  W,  is  used  to  compute  the  independent  components  for  every 
realization,  x(i),  by 


s(i)  =  Wx(i)  (3.8) 

Notice  that  x(i )  is  what  is  observed,  but  s(i)  is  what  needs  to  be  computed. 
Thus,  similar  to  PCA,  the  independent  components  are  determined  indirectly  and  called 
latent  variables  (Hyvarinen,  Oja,  2000).  Moreover,  ICA  is  the  most  common  means  of 
blind  source  separation  (BSS),  where  underlying  factors  arc  determined  “blindly”.  An 
assumption  of  BSS,  which  ICA  necessarily  makes,  is  that  the  number  of  independent 
components  is  limited  to  the  number  of  recorded  signals.  For  instance,  in  the  present  data 
there  were  seven  EEG  and  EOG  signals  recorded,  thus,  ICA  assumes  each  were  a  mixture 
of  seven  independent  components  (Stone,  2002). 

Some  fundamental  assumptions  of  ICA  have  been  mentioned,  but  more 
assumptions  are  necessary  to  make  ICA  feasible.  In  order  to  achieve  components  which 
are  independent,  it  is  necessary  to  assume  the  components  are  non-Gaussian  distributed 
(i.e.  non-normal).  In  this  study,  independent  components  arc  established  by  constructing 
a  mixing  matrix  which  maximizes  non-Gaussianity.  Next,  the  source  signals  are  assumed 
to  propagate  through  a  medium,  in  this  case  brain  tissue,  instantaneously,  before  being 
linearly  mixed  at  the  electrodes.  Finally,  the  source  signals  are  assumed  to  be  stationary. 

ICA  is  made  robust  to  moderate  violations  of  these  assumptions  through  the 
FastICA  algorithm,  detailed  below.  The  major  assumption  which  cannot  be 
compromised,  however,  is  the  number  of  source  signals  assumed  to  exist.  For  the  present 
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purpose,  this  number  is  less  important  than  finding  components  which  are  meaningful  to 
OFS  change  detection  (Makeig  et  al.,  1996;  Vigario  et  al.,  2000). 

Finally,  ICA  differs  fundamentally  from  PCA,  despite  their  similarities.  They  are 
both  multivariate  analyses,  whose  latent  variables  are  linear  combinations  of  observed 
variables.  However,  PCA’s  primary  goal  is  dimensionality  reduction,  whereas  ICA 
creates  a  generative  model  to  reveal  the  underlying  factors  of  a  process.  Additionally, 
PCA  merely  de-correlates  its  latent  variables.  Under  the  assumption  of  normality,  de- 
correlation  is  sufficient  to  achieve  independence.  However,  psychophysiological  signals 
have  consistently  displayed  non-normal  behavior,  implying  that  PCA’s  de-correlation 
does  not  result  in  truly  independent  components.  In  contrast,  ICA  imposes  stricter 
constraints  which  achieve  un-rotatable,  de-correlated,  and  independent  components 
(Stone,  2002). 

Pre-Processing  the  Data  for  Independent  Component  Analysis  Before  conducting 
ICA,  the  data  must  be  pre-processed.  The  data  are  made  to  be  zero  mean,  resulting  in 
independent  components  which  are  also  zero  mean;  this  is  done  to  facilitate  the 
estimation  of  A.  After  it  is  estimated,  the  mean  is  added  back  to  the  independent 
components. 

Next,  the  data  are  whitened  to  make  A  orthogonal  which  significantly  reduces  the 
number  of  parameters  ICA  must  estimate.  Whitening  linearly  transforms  x(£)  into  a  new 
vector,  x(i),  that  is  uncorrelated  and  has  unit  variance.  Whitening  can  be  done  through 
eigenvalue  decomposition  (EVD).  Using  the  same  notation,  x(i)  is  computed  by 

x(i)  =  UL~^2U'x(i)  (3.9) 

This  step  can  also  be  used  to  dimensionally  reduce  the  independent  components,  which 
serves  to  decrease  noise  and  prevent  overfitting  (Hyvarinen,  Oja,  2000). 

The  FastICA  Algorithm  Once  the  data  are  pre-processed,  an  efficient  algorithm  called 
FastICA  is  used  to  estimate  W.  FastICA  determines  the  columns  of  W,  denoted  as  wh 
onc-by-one,  by  maximizing  the  non-Gaussiantity  of  the  projection,  Wj'x(i).  There  arc 
several  measures  of  non-Gaussianity,  among  the  most  common  is  negentropy. 
Negcntropy  is  an  entropic  metric  that  captures  the  “randomness”  of  a  variable.  Because 
Gaussian  variables  are  the  most  random,  they  have  the  highest  entropy  among  all  random 
variables.  Therefore,  FastICA  indirectly  maximizes  non-Gaussianity  by  pursuing  the  least 
entropic  variables  as  measured  by  negentropy.  As  mentioned,  the  non-Gaussianity  of 
components  is  essential  to  achieve  their  independence.  An  approximation  of  negentropy, 
7(y),  is  computed  by 


J(y)  =  (E[G(y)]-E[G(v)])2  (3.10) 

where  G  is  a  non-quadratic  function,  y  and  v  arc  both  zero  mean  and  unit  variance 
variables,  and  v  is  Gaussian  distributed.  In  this  metric,  if  y  is  also  Gaussian  then 
negentropy  is  zero;  any  other  distribution  of  y  results  in  positive  negentropy,  with 
magnitude  proportional  to  y’s  deviation  from  Gaussianity.  In  the  present  project,  the  non- 


16 


quadratic  function  was  specified  as,  G(y)  =  y4,  rendering  (3.10)  a  kurtosis-based 
approximation  of  negentropy. 

The  FastICA  algorithm  is  presented  in  Table  3.2.  Steps  2.3  through  2.5  iteratively 
refine  a  randomly  chosen  vector,  Wj,  until  it  converges.  Convergence  occurs  when  a  new 
Wj  points  in  the  same  direction  as  the  vector  from  the  previous  cycle  of  steps  2.3  and  2.4. 
Once  convergence  occurs,  step  2.8  employs  a  deflation  procedure  to  prevent  different  Wj 
vectors  from  reaching  the  same  maximum;  deflation  de-correlates  the  outputs  of 
wi;t(i),W2*(i),  ...,wn'x(i).  Finally,  step  2.9  renormalizes  Wj  after  deflation.  The 
algorithm  iterates  until  W  is  fully  defined. 

Independent  Component  Analysis  to  Facilitate  Real-Time  OFS  Change  Detection 

In  past  research,  ICA  has  been  conducted  on  physiological  signals  for  two  purposes:  to 
identify  and  extract  artifacts,  such  as  eye  blinks,  and  to  facilitate  the  analysis  of  event- 
related  potentials  (ERP).  There  are  many  artifacts  in  psychophysiological  signals,  such  as 
neck  muscle  activity,  eye  blinks,  heart  rate,  and  line  noise  originating  from  recording 
equipment.  In  some  cases,  the  artifacts’  amplitude  exceeds  the  brain  activity,  the  very 
thing  intended  to  be  measured.  Therefore,  it  is  critical  to  identify  and  extract  these 
artifacts  before  using  psychophysiological  features  to  detect  change  in  OFS.  ICA  can 
automatically  extracts  most  of  the  prominent  artifacts,  producing  less  complex  and  less 
noisy  source  signals  (Vigario  et  al.,  2000). 

ICA  is  also  commonly  used  to  analyze  ERPs  in  EEG  signals.  ERPs  are  spikes  in 
brain  activity  that  result  from  the  onset  of  a  stimulating  event,  like  hearing  a  “beep”.  It  is 
believed  that  ERPs  originate  from  different  sensory  systems  in  the  brain,  corresponding 
to  the  sense  being  stimulated,  like  sight,  sound,  and  touch.  To  identify  the  location  of 
these  sensory  systems,  experiments  are  conducted  which  simultaneously  stimulate 
multiple  senses  of  a  subject  while  recording  their  EEG.  The  ERPs  are  contained  in  the 
EEG,  but  are  mixed  upon  recording.  ICA  is  then  applied  to  “un-mix”  these  signals  and 
generate  independent  components  corresponding  to  the  ERPs  which  originated  from  each 
sensory  system  (Vigario  et  al.,  2000;  Makeig  et  al.,  1996). 

To  the  best  of  our  knowledge,  ICA  has  yet  to  be  employed  to  facilitate  the  online 
change  detection  of  OFS.  Nevertheless,  ICA  is  a  promising  candidate  for  this  application 
since  it  has  been  effective  at  artifact  extraction  and  generating  independent  components 
that  correspond  to  external  stimuli.  The  objective  is  to  use  ICA  to  extract  source  signals 
which  are  free  of  artifacts  and  responsive  to  varying  task  load. 
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Table  3.1  FastICA  algorithm 


1.  Initialize: 

1.1.  7=1 

1 .2.  n  is  the  number  of  independent  components 

2.  Run  Algorithm: 

2.1.  fory<n 


2.2. 

2.3. 


2.4. 


2.5. 

2.6. 

2.7. 

2.8. 


2.9. 


2.10. 


Choose  a  random  vector  Wj 


Wj+  =  £'[A:(i)6,(w;'A:(i))]  —  E[G"(wj'x(i))]wj 
Note:  expectations  are  estimated  by  sample  means. 


w, 


WJ  = 


w, 


if  not  converged,  go  back  to  2.3 
else  if  7=1  go  back  to  2. 1 

end  if/else 

7-i 


wj  =  wj  ~Yj  WJ' 

k=  1 


WkWk 


w, 


w;  = 


Wj  Wj 


7=7+1 


2.11.  end  for 


3.  Output:  the  inverse  mixing  matrix,  W 
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Analyzing  the  Independent  Components  For  each  subject,  all  seven  EEG  and  EOG 
signals  were  submitted  to  ICA  and  the  resulting  independent  components  were  analyzed. 
In  every  case,  at  least  two  independent  components  correlated  directly  with  eye  activity 
measured  at  VEOG  and  HEOG;  this  is  consistent  with  the  visual  demands  of  the  UAV 
tasks.  This  phenomenon  repeated  itself  when  only  EEG  signals  were  submitted  to  the 
ICA,  hence  providing  strong  evidence  that  eye  activity  contaminated  the  EEG  data. 
Figure  3.1  depicts  two  independent  components  extracted  from  the  five  EEG  signals  of 
E01;  it  clearly  displays  their  correlation  to  the  VEOG  and  HEOG  signals  that  were  not 
included  in  the  analysis. 

The  independent  components  that  were  highly  correlated  w  ith  VEOG  and  HEOG 
signals  were  deemed  “artifact”  signals  and  were  excluded  from  further  analysis.  The 
remaining  components  were  empirically  analyzed  to  determine  what  brain  process  they 
represented.  This  was  accomplished  by  constructing  a  frequency  spectrum  for  each 
component,  to  reveal  the  frequencies  that  dominated  them.  Again,  the  DFT,  as  shown  in 
(2.3),  was  used  to  convert  the  components  from  the  time  domain  to  the  frequency 
domain. 

Figure  3.2  depicts  the  frequency  spectrums  of  four  independent  components  from 
A01.  Each  component  was  partitioned  in  HL,  ML,  and  LL  sets  and  submitted  to  the  DFT 
in  three-second  epochs.  The  results  from  the  DFT  were  averaged  across  all  epochs  to 
form  three  frequency  spectrums,  one  for  each  task  load.  As  shown,  each  independent 
component  exhibited  unique  frequency  spectrums,  some  of  which  varied  by  task  load. 
For  instance,  plot  (c)  displays  the  strong  influence  10-12  Hz  waves  had  on  the  third 
independent  component.  Moreover,  the  power  of  these  waves  varied  by  task  load. 

In  fact,  every  subject  had  at  least  one  component  that  w'as  strongly  influenced  by 
8-12  Hz  waves,  corresponding  to  the  alpha  waveband.  Recall  that  this  waveband 
correlates  to  relaxation;  as  a  subject  becomes  more  relaxed,  their  alpha  power  increases. 
This  is  consistent  with  plot  (c)  of  Figure  3.2,  where  alpha  power  is  highest  in  the  LL  and 
lowest  in  the  HL. 

Finding  alpha’s  presence  in  the  psychophysiological  signals  of  cognitive  tasks  is 
not  surprising,  but  revealing  its  prominence  with  such  clarity  is  rare.  Alpha  power  is 
frequently  obscured  by  noise  when  signals  are  mixed  and  recorded  at  the  scalp.  ICA 
segregated  these  noisy  signals  from  the  source  signals  and  clearly  displayed  alpha’s 
contribution.  Makeig  et  al.  confirmed  this  result,  reporting  that  ICA  revealed  “alpha 
activity  (near  10  Hz)  not  obvious  in  the  EEG  data”  (1996).  This  highlights  another  result, 
that  alpha  power  most  often  peaked  around  10  Hz  in  the  independent  components.  This 
particular  frequency’s  significance  has  been  established  in  other  research,  where  it  has 
been  shown  to  correlate  with  performance  on  cognitive  tasks  (Makeig,  Inlow,  1993). 
Research  has  also  shown  10  Hz  to  be  associated  with  spatial  tasks,  more  so  than  verbal 
ones  (Gevins,  1997).  This  is  consistent  with  the  present  results,  since  the  UAV  tasks  were 
largely  spatial  in  nature. 

Analyzing  the  frequency  spectrum  of  a  subject’s  independent  components  may 
prove  useful  in  detecting  changes  of  OFS.  An  abundance  of  research  has  employed  the 
frequency  domain  of  EEG  signals  for  classification,  most  commonly  to  analyze  the 
traditional  wavebands  of  delta,  theta,  alpha,  beta  and  gamma.  Perhaps  a  more  effective 
strategy  would  be  to  use  independent  components  to  derive  subject-specific  and  task- 
specific  wavebands.  An  analysis  similar  to  the  one  depicted  in  Figure  3.2  could  be  done 
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for  each  subject  to  identify  frequency  bands  which  arc  sensitive  to  changes  in  task  load. 
Wilson  and  Fisher  report  using  subject-specific  wavebands  derived  from  latent  variables 
to  increase  classification  accuracy  by  29%  over  traditional  wavebands  (1995).  In  another 
case,  estimating  the  error  rates  on  task  performance  was  improved  significantly  by  using 
a  subject’s  entire  EEG  spectrum  over  predefined,  narrow  wavebands  (Jung  et  al.,  1997). 

However,  one  shortcoming  of  this  method  is  the  subjective  process  of  selecting 
frequency  bands  from  the  independent  components.  It  requires  expert  knowledge  to 
identify  which  components  represent  artifacts,  to  understand  how  brain  activity  responds 
to  different  tasks,  and  finally,  to  decide  which  frequency  bands  are  sensitive  to  task  load. 
These  ambiguities  can  be  mitigated  through  dimensionality  reduction,  where  multiple 
psychophysiological  signals  are  reduced  to  result  in  a  pair  of  independent  components. 
From  this  reduction,  an  independent  component  may  emerge  with  properties  which  arc 
easily  identifiable  and  robust  across  subjects.  The  next  section  discusses  this  approach. 

Dimensionally  Reduced  Independent  Components  During  the  pre-processing  for 
ICA,  eigenvalue  decomposition  (EVD)  was  used  to  whiten  the  data.  EVD  is  equivalent  to 
PCA,  and  like  PCA,  it  can  dimensionally  reduce  the  psychophysiological  signals  to  a 
subset  of  important  components.  Dimensionality  reduction  prevents  overfitting  of  the 
independent  components  and  decreases  noise,  not  to  mention,  reduces  the  vast  number  of 
features  into  a  small,  yet  crucial  subspace.  Most  importantly,  dimensionality  reduction 
reveals  the  independent  components  responsive  to  task  load  with  less  ambiguity. 

EVD  was  used  to  reduce  the  seven  EEG  and  EOG  signals  to  a  pair  of  signals  for 
submission  to  ICA.  Reducing  the  signals  from  seven  to  two  was  justified  by  a  Pareto 
chart  analysis  which,  for  every  subject,  displayed  a  clear  break  at  two  when  plotting  the 
eigenvalues  in  descending  order.  The  resulting  pair  of  signals  was  submitted  to  ICA 
which  generated  two  independent  components.  These  components  were  then  plotted  in 
and  empirically  assessed  for  any  correlation  they  had  w  ith  task  load. 

This  method  yielded  consistent  results  when  conducted  on  each  trial.  In  every 
ease,  one  independent  component  exhibited  clear  variation  with  changing  task  load,  while 
the  other  was  mostly  noise.  Figure  3.3  displays  the  results  of  A01,  where  (a)  is  the 
component  responsive  to  task  load  and  (b)  is  the  noise  component.  Notice  that  the 
independent  component  in  (a)  became  attenuated  during  the  HL  and  ML,  compared  to  its 
variation  in  the  LL.  The  second  component  in  (b)  exhibited  similar  behavior,  but  was  not 
consistent  for  every  task  and  was  obscured  by  noise.  The  utility  of  the  first  component  for 
real-time  OFS  change  detection  is  promising  since:  one  exists  for  each  subject,  its 
behavior  is  distinguishable  by  task  load,  and  finally,  the  ICA  method  is  automatic  and 
leaves  little  error-prone  subjectivity. 

Further  investigation  of  the  dimensionally  reduced  independent  components 
disclosed  strong  correlations  with  the  VEOG  and  HEOG  signals.  Previously,  the 
independent  components  correlated  with  VEOG  and  HEOG  were  excluded  from  further 
analysis  since  they  were  deemed  “artifacts”.  However,  in  the  present  analysis  these 
components  were  extracted  as  the  source  signals  for  all  the  EEG\EOG  data,  thus  they 
were  retained. 
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Figure  3.1  (a)  IC1  (b)  VEOG  (c)  IC2  and  (d)  HEOG  signals  of  E01 
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Figure  3.2  Frequency  spectrum  of  independent  components  from  A01,  (a)  IC1  (b)  IC2  (c)  IC3  (d)  IC4 
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Table  3.3  and  Table  3.4  contain  the  correlations  between  the  original 
psychophysiological  signals  and  the  first  and  second  independent  components, 
respectively.  The  correlations  are  averaged  across  trails  by  subject.  Subject  F  is 
considered  separately  since,  when  his  components  were  plotted,  they  exhibited  unique 
behavior.  Subject  F  has  consistently  been  the  most  difficult  to  accommodate  in  a  variety 
of  methods,  and  the  correlation  analysis  provided  clues  as  to  why. 


Table  3.3  Correlation  analysis  of  dimensionally  reduced 
first  independent  component 


VEOG 

HEOG 

Fz 

f7 

Pz 

t5 

02 

Subjects  A,  E 

0.903 

0.107 

0.835 

0.385 

0.389 

0.297 

0.179 

Subject  F 

0.977 

0.075 

-0.380 

-0.137 

-0.261 

-0.228 

-0.174 

Table  3.4  Correlation  analysis  of  dimensionally  reduced 
second  independent  component 


— 

VEOG 

HEOG 

Fz 

f7 

Pz 

t5 

o2 

Subjects  A,  E 

0.373 

0.992 

0.067 

0.860 

0.281 

0.036 

0.475 

Subject  F 

0.213 

-0.991 

-0.104 

-0.870 

0.045 

-0.159 

0.227 

In  Table  3.3,  notice  the  strong  correlation  the  first  independent  component  had 
with  the  VEOG  signal  across  subjects.  Another  notable  correlation  for  subjects  A  and  E, 
is  with  the  EEG  signal  recorded  at  Fz.  However,  this  is  not  true  for  subject  F,  who 
averaged  less  than  half  the  correlation  at  Fz  than  subjects  A  and  E. 

The  second  independent  component,  in  Table  3.4,  correlated  with  one  EOG  and 
one  EEG  signal,  similar  to  the  first  component.  But  in  this  case,  the  correlation  was  with 
HEOG  and  F7,  respectively.  Notice  that  subject  F’s  behavior  was  not  appreciably 
different  than  the  other  subjects. 

This  analysis  elucidates  several  important  findings,  principally,  that  the  source 
signals  of  the  EEG\EOG  were  strongly  correlated  to  eye  activity.  This  is  useful 
knowledge,  as  otherwise,  VEOG  and  HEOG  contaminated  signals  may  be  discarded  as 
artifacts  rather  than  treated  as  important.  As  mentioned,  one  of  the  two  independent 
components  always  varied  with  task  load  and  the  other  was  mostly  noise.  For  subjects  A 
and  E,  the  independent  component  which  best  varied  with  task  load  was  the  first.  The 
correlations  in  Table  3.3  indicate  that  this  component  was  a  mixture  of  signals  from 
VEOG  and  Fz.  Recall  that  the  Fz  electrode  mostly  records  theta  waves  which  arc  sensitive 
to  changes  in  cognitive  load.  This  was  true  except  for  subject  F,  who  did  not  have  a 
strong  Fz  presence  in  the  first  component.  In  fact,  subject  F’s  first  component  was  noisy 
and  it  was  his  second  component,  the  one  dominated  by  HEOG  and  F7,  that  best  varied 
with  task  load.  Two  things  can  be  inferred  from  this:  first,  ICA  can  reveal  the  source  of 
between-subject  differences,  and  second,  those  differences  can  be  overcome  with  at  least 
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one  independent  component  that  is  sensitive  to  task  load.  Henceforth,  the  independent 
component  that  best  varied  with  task  load,  post-dimensionality  reduction,  is  referred  to  as 
the  TVIC  (task-varying  independent  component). 


4.  Real-Time  Metrics  for  Psychophysiological  Signals  and  Algorithms 
for  Temporal  Change  Detection  in  Time  Series  Data 

The  Peak  Detection  Method  Now  that  a  robust  signal,  the  TVIC,  can  be  identified  for 
each  subject,  a  quantitative  metric  is  necessary  to  characterize  the  signal  for  each  epoch; 
this  metric  will  be  used  to  facilitate  the  online  change  detection  of  OFS.  The  TVIC 
plotted  in  (a)  of  Figure  3.3  suggests  that  a  metric  of  variance  might  best  characterize  the 
signal’s  response  to  task  load,  since  variance  noticeably  increased  during  the  LL  and 
drastically  reduced  during  the  ML  and  HL.  Ideally,  monitoring  a  metric  like  variance, 
computed  for  each  epoch,  would  result  in  a  signal  with  less  noise  and  better 
discrimination  between  task  loads  than  the  raw  TVIC.  In  essence,  the  magnitude  of  the 
metric  would  “peak”  when  the  tasks  are  incurred,  thus  signifying  changes  to  higher 
cognitive  loads. 

Although  variance  is  an  obvious  choice,  there  arc  other  quantitative  metrics  to 
characterize  a  signal  in  real-time,  such  as  the  third  moment  (skewness),  the  fourth 
moment  (kurtosis),  and  finally,  entropy  based  metrics.  These  metrics  may  also  prove 
sensitive  to  the  onset  of  tasks  and  even  result  in  superior  change  detection. 

An  objective  method  was  needed  to  evaluate  and  compare  each  metric’s  ability  to 
produce  consistent  task-induced  peaks.  One  such  technique  is  called  the  peak  detection 
method  (PDM).  The  PDM  compares  a  measure  Xi  at  epoch  /  with  two  thresholds,  an 
absolute  threshold  defined  by  parameter  y  and  a  relative  threshold  defined  by  parameter 
5;  the  comparisons  are  made  to  determine  if  a  peak  has  occurred. 

For  the  PDM,  each  metric  was  computed  by  a  sliding  window  so  that 
irregularities  w'ere  smoothed  and  peaks  were  induced  during  the  tasks.  The  PDM  was 
used  to  objectively  evaluate  which  metric  most  consistently  peaked  for  all  8  HL  and  ML 
tasks,  while  generating  the  fewest  false  alarms.  Hence,  the  results  of  the  PDM  indicated 
which  metric  best  detected  change  in  OFS.  The  algorithm  for  the  PDM  is  presented  in 
Table  4.1  (Yu,  2009). 
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Figure  3.3  Dimensionally  reduced  independent  components  of  A01,  (a)  IC1  (b)  IC2 
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Review  of  Metrics  for  Characterizing  Psychophysiological  Signals  The  simplest 
metric  for  characterizing  the  TV1C  is  the  standard  deviation.  The  standard  deviation  is  the 
positive  square  root  of  the  signal’s  second  moment,  defined  as 


where  Mj  is 


a  =  JW2 

£”=i(*i  -/0; 


(4.1) 


(4.2) 


H  is  the  mean  of  the  signal,  generally  estimated  by  the  sample  mean,  x,  computed  from  a 
sample  of  size  n.  For  the  standard  deviation  calculation,  n  is  traditionally  reduced  by  one 
in  the  denominator  of  (4.2)  (Montgomery,  2009). 

A  related  metric  is  the  third  moment  of  the  TV1C,  termed  skewness  when 
normalized,  which  measures  the  asymmetric  nature  of  a  distribution.  The  TVIC  in  plot 
(a)  of  Figure  3.3  displays  clear  skewness  in  the  distribution  of  its  amplitudes  during  the 
LL,  but  much  less  skewness  during  the  tasks.  Skewness  is  defined  as 

M3 

Y  Ml12  (4.3) 


(Montgomery,  2009). 

The  fourth  moment  may  also  prove  a  valuable  metric  as  it  characterizes  the  shape 
of  a  distribution,  discriminating  distributions  that  are  tall  and  skinny  from  those  that  are 
short  and  stout.  When  standardized,  the  fourth  moment  is  termed  kurtosis,  and  because  it 
is  fourth  order,  it  is  always  positive.  Kurtosis  is  defined  as 


k 


(4.4) 


(Montgomery,  2009). 

The  next  two  metrics  are  both  measures  of  entropy,  where  entropy  quantifies  the 
“randomness”  of  the  TVIC.  The  more  random  the  data,  the  higher  entropy  it  will  have, 
conversely,  the  more  ordered  the  data,  the  less  entropy  it  will  have;  data  in  perfect  order 
has  zero  entropy.  In  recent  research,  entropy  measures  of  psychophysiological  data  have 
been  found  to  correlate  with  changing  states  of  vigilance,  in  particular,  entropy  increased 
with  increasing  vigilance  (Bruzzo  et  al.,  2008;  Zhang  et  ah,  2009).  This  result  suggests 
entropy  metrics  may  be  able  to  detect  changes  in  OFS. 

The  first  entropy  metric  considered,  called  sample  entropy  (SampEn),  is  a  non¬ 
linear  metric  that  measures  the  regularity  of  a  time  series  (Richman  et  ah,  2004).  SampEn 
was  chosen  because  it  is  computationally  efficient  and  can  be  applied  to  short,  noisy  data 
typical  of  the  TVIC.  In  addition,  SampEn  is  an  improvement  over  the  more  traditional 
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time  series  entropy,  approximate  entropy  (ApEn);  SampEn  exhibits  less  bias  towards 
lower  entropy,  improves  the  measure’s  relative  consistency  between  datasets,  and  is  more 
robust  with  regards  to  varying  record  lengths  (Richman,  Moorman,  2000). 

SampEn  measures  the  regularity  of  a  time  series  by  computing  the  conditional 
probability  that  two  arbitrarily  similar  epochs  of  size  m,  remain  arbitrarily  similar  for  the 
next  point  in  the  series.  “Arbitrary”  is  used  here  to  emphasize  a  degree  of  error  allowed 
when  classifying  two  epochs  as  similar.  This  is  controlled  by  the  parameter  r,  typically  a 
factor  between  .1  and  .25  of  the  record’s  standard  deviation.  By  SampEn’s  definition,  a 
perfectly  ordered  series  will  have  a  conditional  probability  of  one,  corresponding  to 
entropy  of  zero.  Time  series  with  less  order  will  have  positive  entropies  and  conditional 
probabilities  less  than  one  (Richman  et  al.,  2004).  The  algorithm  to  compute  SampEn  is 
in  Table  4.2  (Alcaraz,  Rieta,  2008). 

The  second  entropy  metric,  termed  Kullback  Liebler  Divergence  (KLD), 
measures  the  relative  entropy  between  two  probability  distributions,  p  =  { pk }  and 
q  =  {qk}.  In  the  present  study,  frequency  distributions  were  analyzed  instead  of 
probability  distributions,  where  A:  is  a  frequency  and  pk  is  the  normalized  density  at  k, 
similarly  for  qk  (Quiroga  et  al.,  2000).  The  KLD  computed  the  relative  entropy  between  a 
static  frequency  distribution  constructed  from  the  LL  data,  and  a  frequency  distribution 
that  changed  with  a  sliding  window  over  the  TVIC.  The  hypothesis  was  that  the 
frequency  distribution  from  the  HL  and  ML  differed  markedly  from  the  LL  distribution, 
thus  KLD  would  peak  during  these  tasks. 

KLD  is  always  positive  and  its  magnitude  is  proportional  to  the  difference 
between  distributions;  it  is  zero  for  identical  distributions.  KLD  is  computed  by 


(4.5) 


Results  of  the  PDM  on  Various  Metrics  Several  different  metrics  were  presented, 
each  with  the  ability  to  characterize  the  task-varying  nature  of  the  TVIC.  The  PDM  was 
used  to  objectively  evaluate  which  of  the  metrics  most  consistently  peaked  for  tasks.  This 
section  presents  those  results,  and  identifies  the  metric  which  is  best  suited  for  online 
OFS  change  detection. 

The  standard  deviation,  skewness,  and  kurtosis  metrics  were  all  computed  by  a 
sliding  window  over  each  subject’s  TVIC.  The  span  of  the  window  was  18  seconds, 
providing  an  update  of  the  metric  each  second.  The  length  of  the  span  was  chosen  for 
optimal  peak  production,  as  18  seconds  corresponds  to  the  duration  of  the  tasks;  hence,  at 
one  point  the  span  of  the  window  would  precisely  overlap  an  entire  task,  resulting  in  the 
most  extreme  realization  of  the  metric,  i.e.,  the  largest  peak.  The  span  may  seem 
extensive,  but  it  was  necessary  to  smooth  the  metrics  and  resolve  any  task-induced 
trends. 

For  the  following  analyses,  the  S  and  y  parameters  of  the  PDM  were  optimized 
with  respect  to  each  subject  and  metric  under  evaluation.  This  strategy  is  consistent  with 
operational  requirements,  as  between-subjeets  variation  is  too  profound  to  generalize 
parameter  settings.  An  example  of  the  standard  deviation  of  the  TVIC  from  E02  is  plotted 
in  Figure  4.1;  the  standard  deviation  is  negated  for  interpretability  and  the  peaks  detected 
via  the  PDM  are  indicated  by  cross-hairs.  All  8  tasks  were  detected. 
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Table  4.1  PDM  algorithm 


1. 

Initialize: 

1.1. 

Testing  set,  5  =  { x1,x2 , ...  ,xn } 

1.2. 

8,  y  are  user  specified 

1.3. 

Lmax  ~  00  5  I max—  0 

1.4. 

^min  ~  00  5  I min—  1 

1.5. 

B  =  1 

1.6. 

P=0,I=0 

1.7. 

i=  1 

2. 

Run  Algorithm: 

2.1. 

while  /  <  n 

2.2. 

if  xi  >  Lmax 

2.3. 

L max  ~  xi ,  I max  ~  * 

2.4. 

end  if 

2.5. 

if  Xi  <  Lmin 

2.6. 

Lmin  ~  xi 

2.7. 

end  if 

2.8. 

if  B  =1 

2.9. 

if  xi  <  Lmax  8  AND  Lrnax  ^  Y 

2.10. 

P^PuLmax 

2.11. 

I  <  /  u  lmax 

2.12. 

B=  0 

2.13. 

Lmin  ~  xi 

2.14. 

end  if 

2.15. 

else  if  Xi  >  Lmin  +  8 

2.16. 

Lmax  ~  xi 

2.17. 

Imax~i 

2.18. 

B  =1 

2.19. 

end  if/else 

2.20. 

/  =  /  +1 

2.21. 

end  while 

3. 

Output:  peaks  and  their  indices  are  contained  in  sets  P  and  /,  respectively 
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Table  4.2  SampEn  algorithm 


1. 

Initialize: 

1.1. 

N  is  the  length  of  the  record 

1.2. 

m  and  r  are  user  defined  parameters 

1.3. 

Form  vectors  of  length  m,  xm{l),  xm(2 ), ... 

,xm(N  -  m  +  1) 

1.3.1.  xm(i )  =  [x(i),x(i  +  1),  ...,x(i  +  m  - 

- 1)]  v  / 

1.4. 

Bh  A[,  Bm,  Am,  D  are  initialized  to  zero 

1.5. 

i=  1  ,j  =  1 

2. 

Run  Algorithm: 

2.1. 

while  i  <  N  —  m  +  1 

2.2. 

b  =  0 

2.3. 

while  j<N  —  m+l  AND  j  ^  i 

2.4. 

Select  two  vectors  xm(i)  and  xm(J ) 

2.5. 

D  =  max  1  x(i  +  k)  — 

fc=0,..,m-l 

x(j  +  k)  | 

2.6. 

if  D  <  r 

2.7. 

b  =b  +  1 

2.8. 

end  if 

2.9. 

end  inner  while 

2.10. 

1 

D  — 

U 

1  N  —  m  —  1 

2.11. 

end  outer  while 

2.12. 

1 

pm  _ 

N-m 

.  V  p. 

N  —  m  Z_i 

i= 1 

2.13. 

Repeat  steps  1-3,  except  for  vectors  of  length  m  +  l.j  ranges  from  1  to 

N  —m  and  replace  b  with  a,  and  B  with  A 

3. 

Output: 

SampEn  =  —In 

r/Pi 

Bm 
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SampEn  was  also  computed  by  a  sliding  window,  in  this  case,  with  a  span  of  10 
seconds.  The  span  was  shorter  because  a  second  window  was  applied  in  addition  to  the 
first.  SampEn  exhibited  brief,  erratic  outliers  whose  false  peaks  were  not  sufficiently 
smoothed  by  a  single  window.  To  compensate,  a  second  window  was  applied  over  the 
first-windowed  TVIC,  and  for  each  set  of  10  points,  it  computed  the  area  under  the  curve, 
i.e.,  approximated  the  integral.  The  double-windowing  smoothed  false  peaks,  as  their 
short  duration  was  not  associated  with  large  area  under  the  curve.  Figure  4.2  displays  an 
example  of  this  for  A02.  In  (a),  the  first  sliding  window  of  SampEn  was  applied  to  the 
TVIC.  Notice  the  presence  of  several  short,  narrow  outlier  peaks  in  contrast  to  the  wider 
peaks  exhibited  during  each  task.  Plot  (b)  is  the  integral  from  the  second  sliding  window 
over  the  first  windowed  signal  in  (a);  it  demonstrates  the  dampening  of  outlier  peaks  and 
the  correct  detection  of  all  8  tasks  with  no  false  alarms.  The  parameters  from  SampEn 
were  set  as  follows:  r  to  .2  times  the  standard  deviation  of  the  record,  m  to  2,  and  the 
record  length,  N,  to  2000  (Richman  et  al.,  2004). 

The  last  metric  analyzed  was  the  KLD  which  computed  the  relative  entropy 
between  two  frequency  distributions.  The  static  frequency  distribution  was  derived  by 
averaging  the  frequencies  of  each  10  second  epoch  from  the  LL  of  an  entire  trial’s  TVIC; 
as  before,  the  frequencies  were  computed  via  the  DFT.  The  second  distribution  was 
computed  in  real-time  from  a  sliding  window  over  the  subject’s  TVIC.  The  density  of 
each  frequency  distribution  was  normalized,  so  the  area  under  the  distribution  was  one. 
After  normalization,  both  the  static  baseline  distribution  and  the  real-time  distribution 
were  submitted  to  the  KLD,  for  every  second,  to  form  a  relative  entropy  signal. 

The  span  of  the  sliding  window  for  the  KLD  was  again  set  to  10  seconds,  because 
as  with  SampEn,  the  KLD  necessitated  a  second  window  to  dampen  outlier  peaks.  As 
before,  the  second  window  slid  over  the  initial  windowed  signal  and  computed  the  area 
under  the  curve  for  each  set  of  10  points.  Figure  4.3  is  an  example  for  E02  where  plot  (a) 
is  the  KLD  windowed  over  the  TVIC  and  plot  (b)  is  the  integral  from  the  second  sliding 
window  over  the  first  windowed  signal  in  (a).  Similar  to  Figure  4.2,  the  double 
windowing  smoothed  outliers  while  amplifying  the  task-induced  peaks  to  detect  7  of  the 
8  tasks. 

The  overall  results  of  the  PDM  analysis  are  in  Table  4.3.  The  proportion  of 
correctly  detected  tasks  is  listed  for  each  metric  by  trial.  Average  accuracy  and  false 
alarm  rate  is  also  reported  in  bold  for  each  metric.  According  to  the  results,  the  standard 
deviation  was  the  most  proficient  at  producing  peaks  for  each  task.  Using  the  standard 
deviation,  the  PDM  detected  93.8%  of  the  tasks  with  only  .5  false  alarms  per  trial  on 
average.  The  standard  deviation  dominated  all  other  metrics,  i.e.,  for  no  subject  or  trial 
did  another  metric  perform  better.  SampEn  also  had  strong  performance,  detecting  81.3% 
of  tasks  on  average,  but  it  failed  for  subject  F.  The  next  most  consistent  metric  was 
kurtosis,  with  68.8%  of  the  tasks  correctly  detected  on  average.  Recall  that  kurtosis 
quantifies  the  shape  of  a  distribution,  meaning  that  with  some  success,  task  loads  can  be 
distinguished  by  the  shape  of  the  TVIC’s  amplitude  distribution.  KLD  also  utilized 
distributions,  in  this  case,  to  quantify  differences  in  entropy.  The  KLD  metric  detected 
64.6%  of  tasks  on  average  in  the  PDM. 
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Table  4.3  Results  of  the  PDM  on  various  metrics 


Std. 

Deviation 

Skewness 

Kurtosis 

SampEn 

KLD 

A01 

1.000 

1.000 

1.000 

1.000 

1.000 

A02 

1.000 

0.250 

0.750 

1.000 

0.750 

E01 

1.000 

0.375 

0.625 

1.000 

0.750 

E02 

1.000 

0.750 

0.625 

1.000 

0.875 

F01 

0.875 

0.625 

0.750 

0.500 

0.250 

F02 

0.750 

0.625 

0.375 

0.375 

0.250 

Average 

0.938 

0.604 

0.688 

0.813 

0.646 

FA/Trial 

0.500 

0.000 

1.250 

0.750 

0.750 

Discussion  of  Metrics  and  the  PDM  The  PDM  was  used  as  a  tool  to  objectively 
evaluate  how  well  each  metric  could  detect  change  in  task  loads.  The  PDM  analysis 
yielded  that  the  windowed  standard  deviation  of  the  TVIC  was  the  best.  By  monitoring 
this  signal  and  metric  for  each  subject,  a  change  detection  scheme  could  accurately 
identify  when  the  subject  changes  state. 

The  entropy  metrics  were  also  competitive  in  detecting  tasks,  but  only  for  subjects 
A  and  E.  Once  again,  subject  F's  psychophysiological  signals  proved  difficult  to 
accommodate,  as  they  did  not  exhibit  a  recognizable  entropy  pattern.  With  data  on  more 
subjects,  the  entropy  metrics  may  prove  excellent  in  characterizing  task-varying  signals, 
on  the  other  hand,  it  may  be  reaffirmed  that  entropies  are  largely  subject-specific. 

Yu  first  introduced  the  PDM’s  application  to  the  UAV  dataset.  He  presented  task 
detection  results  from  the  PDM  on  the  windowed  standard  deviation  of  a  subject’s  F7 
signal  (2009).  The  span  of  the  sliding  window  was  set  to  18  seconds  and  the  PDM 
parameters,  8  and  y,  were  fixed.  In  order  to  make  a  fair  comparison  with  these  findings, 
the  PDM  analysis  was  repeated,  this  time  with  fixed  parameters  on  the  standard  deviation 
of  a  subject’s  TVIC.  The  outcome  of  the  comparison  is  in  Table  4.4. 

The  average  accuracy  using  a  subject’s  TVIC  was  2%  higher  than  using  their  F? 
signal.  While  this  is  not  a  significant  improvement,  the  PDM  on  the  TVIC  yielded  less 
than  half  the  false  alarms  on  average  when  compared  to  the  Fz  signal. 

The  results  support  the  notion  of  using  independent  components  instead  of  the 
raw  psychophysiological  signals  for  OFS  change  detection.  Independent  components  arc 
less  complex  and  less  noisy,  and  thus  result  in  fewer  false  alarms.  Not  only  do  these 
results  confirm  ICA’s  ability  to  extract  artifacts,  they  directly  demonstrate  ICA’s  ability 
to  facilitate  the  detection  of  change  in  a  subject’s  cognitive  state. 

The  PDM  is  a  valuable  tool,  but  its  role  should  not  be  extrapolated  to  that  of  an 
online  change  detector.  In  the  context  of  the  present  study,  a  peak  is  only  defined  once  a 
metric  reverses  direction,  and  begins  descending  back  to  the  LL  baseline.  In  other  words, 
a  peak  is  only  detected  once  the  cognitive  load  returns  to  normal,  long  after  a  task  has 
ended.  Therefore,  the  PDM  cannot  be  a  true  online  change  detector,  as  it  can  only  detect 
tasks  in  retrospect. 
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Figure  4.1  PDM  on  the  standard  deviation  of  the  TVIC  of  E02 
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Figure  4.2  PDM  on  the  SampEn  of  the  TVIC  of  A02  (a)  first  sliding  window  (b)  second  window  of  integral  on  (a) 
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Figure  4.3  PDM  on  the  KLD  of  the  TVIC  of  E02  (a)  first  sliding  window  (b)  second  window  of  integral  on  (a) 
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Table  4.4  Results  of  the  PDM  on 
the  standard  deviation  of  Fz  and  TVIC 


A01 

Fz  TVIC 

1.000  1.000 

0.875  0.750 

1.000  1.000 

0.750  0.875 

1.000  0.875 

0.500  0.750 

A02 

E01 

E02 

F01 

F02 

Average 

0.854  0.875 

2.000  0.833 

FA/Trial 

The  PDM  does  provide  an  interesting  model,  one  that  could  be  modified  into  a 
real-time  scheme  that  detects  a  metric’s  task-induced  trends.  After  identifying  the 
standard  deviation  as  the  optimal  metric  to  characterize  a  subject’s  TVIC,  the  next  step  is 
to  employ  the  metric  in  a  real-time  OFS  change  detector.  The  following  section 
introduces  a  novel  algorithm,  loosely  based  on  the  PDM,  which  monitors  and  detects 
task-induced  trends  instead  of  peaks. 

The  Trend  Detection  Method  A  true  online  classifier  cannot  benefit  from  hindsight, 
meaning  the  information  used  to  classify  an  epoch  of  data  must  come  from  that  epoch 
and/or  the  information  that  preceded  it.  The  PDM  does  not  qualify  as  a  true  online  change 
detector  because  it  requires  information  about  what  transpires  after  an  epoch  in  order  to 
retroactively  detect  peaks.  An  online  classifier  must  also  be  instantaneous,  and  most 
importantly,  it  must  be  accurate. 

The  trend  detection  method  (TDM)  introduced  in  this  section  was  developed  to 
meet  these  requirements  and  to  address  the  deficiencies  of  the  PDM.  The  TDM  assumes 
that  the  ML  and  HL  will  cause  a  windowed  metric  to  monotonically  increase  (or 
decrease)  from  the  LL  baseline.  This  behavior  was  observed  when  monitoring  the 
standard  deviation  of  a  subject’s  TVIC.  The  TDM  identifies  and  detects  these  task- 
induced  trends  in  order  to  detect  real-time  changes  in  OFS. 

The  keys  to  trend  detection  are  to  identify  significant  trends  through  noise  and  to 
do  so  for  non-stationary  signals,  where  trends  can  begin  at  different  magnitudes  over 
time.  The  TDM  addresses  these  challenges  via  an  algorithm,  presented  in  Table  4.,  which 
has  two  main  components:  an  adaptive  threshold  and  a  trend  detector.  The  adaptive 
threshold  adjusts  to  non-stationary  behavior  and  the  second  component,  the  trend 
detector,  identifies  task-induced  trends.  Detecting  trends  accomplishes  two  things:  first, 
coupled  with  information  on  the  signal’s  magnitude,  a  trend  helps  discriminate  between 
true  task-induced  increases  and  irrelevant  transient  spikes.  Second,  depending  on  the 
length  of  the  trend,  tasks  occurring  below  the  threshold  can  be  identified  and  adapt  the 
threshold  accordingly. 

The  TDM  judges  a  signal,  xt,  at  time  i  against  two  criteria:  its  magnitude  relative 
to  the  adaptive  threshold  and  how  well  the  previous  points,  xt _1,xi_2,  ....x^,  have 
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trended.  Trends  are  identified  by  counting  the  number  of  positive  “slopes”  between 
consecutive  points  over  a  period  of  time,  where  r  is  the  count.  A  positive  slope  is  defined 
when  Xi  >  xi_1.  The  slopes  need  not  be  consecutive  since  noise  may  cause  temporary 
reversals  in  the  trend;  reversals  will  only  dissolve  a  trend  if,  after  a  period  of  time,  no 
point  falls  above  the  trend’s  last  point,  say,  Xj.  This  time  period  is  defined  by  the 
parameter  im.  If,  however,  Rum  is  not  exceeded  and  some  point,  xi+n  >  xt-,  then 
r  =  r  +  1  and  the  trend  detection  continues.  When  >  Ttow,  the  second  phase  of  the 
algorithm  is  commenced. 

Once  an  emerging  trend  is  identified,  the  last  value  in  the  trend,  X;,  is  compared  to 
the  threshold,  Z,  to  determine  if  a  task-induced  cognitive  load  has  occurred.  If  x*  >  Z,  the 
trend  is  classified  as  task-induced.  In  contrast,  if  x(-  <  Z,  the  emerging  trend  remains 
unclassified  but  continues  to  be  monitored  until  r  >  Tup ,  where  Tup  is  a  newly  defined 
limit.  In  summary,  the  second  phase  of  the  algorithm  classifies  trends  in  one  of  two  ways: 
the  first  is  when  the  trend  outright  breaches  Z,  and  the  second,  is  when  the  trend  is 
occurring  below  Z,  but  is  sustained  long  enough  to  be  deemed  task-induced.  Tup 
quantifies  “long  enough”  by  accounting  for  the  distance  xt  falls  below  Z;  the  further 
below  Z,  the  longer  trending  must  continue  in  order  to  be  classified.  The  motivation  here 
is  to  account  for  non-stationarity  while  preventing  the  classification  of  trends  not 
associated  with  a  task  load. 

If  a  task-induced  trend  is  classified  below  the  current  Z,  the  threshold  updates  by 
an  EWMA  that  weights  the  magnitudes  of  previously  classified  trends  with  the  most 
recent  trend,  through  a  smoothing  parameter  X.  This  mechanism  adapts  the  threshold  so  it 
can  detect  task-induced  trends  occurring  at  different  magnitudes  in  the  future.  These 
mechanisms  of  the  TDM  are  all  detailed  in  Table  4.5. 

Results  of  the  TDM  on  the  Standard  Deviation  of  TVIC  The  TDM  was  conducted 
on  the  standard  deviation  of  an  18  second  sliding  window  over  each  subject’s  TVIC;  this 
is  the  same  setup  as  the  PDM  analysis.  Normally,  the  parameters  of  the  TDM  are 
calibrated  during  training.  However,  for  the  purpose  of  comparing  results  to  the  PDM,  the 
parameters  were  calibrated  on  the  testing  data.  In  general,  the  parameters  were  set  within 
the  following  ranges:  Tiow  [4,6],  L  [1.2, 1.4],  /?  [8,9],  The  a  and  Rnm  parameters  were 
both  held  constant  at  .3  and  5,  respectively.  Figure  4.4  is  a  plot  of  the  TDM  of  E01;  8 
task-induced  trends  were  detected,  with  one  false  alarm.  Notice  the  instances  a  trend  was 
identified  more  than  once  for  the  same  task;  this  only  indicates  the  TDM  detected  the  task 
early  and  continued  to  detect  further  trending  as  the  task  proceeded. 

The  results  of  the  TDM  analysis  for  each  trial  are  in  Table  4.5.  The  accuracy  is 
defined  as  the  proportion  of  tasks  that  were  correctly  detected.  The  average  accuracy  of 
the  TDM  was  85.4%,  with  over  one  false  alarm  per  trial  on  average.  Lastly,  the  tasks 
were  detected,  on  average,  13  seconds  after  their  onset,  a  comparatively  quick  detection 
when  compared  to  the  half  minute  delays  in  control  charting. 
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Table  4.5  TDM  algorithm 


1.  Initialize: 

1.1.  Testing  set,  S  =  {x1,x2,  ...,xn} 

1.2.  Tlow,  Rum,  A,  /?,  and  L  are  all  defined  during  training 

1.3.  r  =  0,  v=  0 

1.4.  Tup  =  oo+ 

1 .5.  Lmax=  where  /j.x  is  the  mean  of  S  or  a  training  set  X 

1 .6.  Z0  =  Hx  +  Lcrx,  where  ox  is  the  standard  deviation  S  or  of  a  training  set  X 

1.7.  V  =  0,1  =  0 

1.8.  /=  1 

2.  Run  Algorithm: 

2.1.  do 

2.2.  /=/+ 1 

2.3.  \{xi>xi_1 

2.4.  if  Xi  >  Lmax 

2.5.  Lmax-  Xi 

2.6.  r  =  r+l 

2.7.  v  =  0 

2.8.  else  v  =  v  + 1  end  if/else 

2.9.  else  v  =  v  +  1  end  if/else 

2. 1 .  while  (  v  <  Rum  AND  r  <  Tiow  AND  /  <  t ) 

2.2.  if  v>RUm 

2.3.  r  =0 

2.4.  3,max  —  X( 

2.5.  else  if  xt  >  Z  OR  r  >  Tup 

2.6.  V  —  V  U  Xj 

2.7.  I  <—  I  U  / 

2.8.  r  =  0,  v  =  0 

2.9.  Tup  -  00 + 

2.10.  Lmax=fix 

2.11.  Z  =min{  Z0,  Z  (\-A.)  +  ALmax} 
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Tabic  4.5  Results  of  the  TDM  analysis 


Trial 

Proportion  Detected 

A01 

1.00 

A02 

.750 

E01 

1.00 

E02 

.875 

F01 

.750 

F02 

.750 

Average 

.854 

FA/Trial 

1.166 

Discussion  of  the  TDM  Overall,  the  results  indicate  that  the  TDM  detected  changes  in 
cognitive  load  accurately,  quickly,  and  most  importantly,  in  real-time.  Pattern  recognition 
methods  that  classify  epochs  into  discrete  categories  of  eognition  arc  not  eondueive  to 
real-time  change  detection;  it  is  unclear  when  a  subject  changes  state.  In  contrast,  the 
TDM,  in  combination  with  a  subject’s  TVIC,  can  clearly  indicate  when  a  change  has 
occurred. 

The  principle  shortcomings  of  the  TDM  arc  its  inability  to  signal  when  a  task  has 
ended  and  the  number  of  parameters  that  must  be  calibrated.  In  adaptive  aiding,  it  is 
necessary  for  the  TDM  to  not  only  identify  the  beginning  of  a  task,  but  also  indicate  the 
task’s  conclusion  so  aiding  ean  be  withdrawn.  This  capability  ean  be  added  to  the 
algorithm  by  monitoring  trends  that  reverse  to  the  baseline  onee  a  task  completes.  The 
second  drawback  is  the  number  of  parameters  that  need  to  be  set  for  the  TDM,  five  in  all. 
In  the  present  analysis,  only  three  parameters  required  adjustment  across  trials,  but  even 
setting  these  parameters  is  not  straightforward.  This  process  could  be  formalized  by 
conducting  a  design  of  experiments  on  training  data,  scientifically  determining  the 
optimal  parameter  settings  for  each  subject. 

These  suggestions  to  remedy  the  shortcomings  of  the  TDM  wrere  not  pursued 
further,  since  the  scope  of  this  study  is  a  proof-of-eoneept  rather  than  an  optimization  of 
any  single  method.  The  goal  of  employing  adaptive  aiding  from  psychophysiological 
signals  will  require  more  than  a  single  method.  The  TDM  and  the  TVIC  have  been 
developed  to  augment  the  arsenal  of  sehemes  which  identify  changing  cognitive  load,  and 
it  is  hoped  that  in  concert  with  these  established  methods,  future  adaptive  aiding  schemes 
can  be  realized. 
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Figure  4.4  TDM  on  the  standard  deviation  of  the  TVIC  of  E01 


APPENDIX  A 

TASK  LOAD  MEANS  ANALYSIS 

Table  A.  1  Means  of  features  across  task  loads  for  A01 


Theta 

Alpha 

Beta 

LL 

1.182 

0.637 

0.988 

ML 

1.246 

0.609 

0.809 

HL 

1.213 

0.572 

0.614 

Table  A. 2  Means  of  features  across  task  loads  for  E01 


Theta 

Alpha 

Beta 

LL 

1.685 

0.569 

1.023 

ML 

1.564 

0.447 

0.654 

HL 

1.493 

0.635 

0.870 

Table  A. 3  Means  of  features  across  task  loads  for  F01 


Theta 

Alpha 

Beta 

LL 

0.987 

0.332 

0.585 

ML 

0.925 

0.277 

0.414 

HL 

1.026 

0.263 

0.439 

Table  A.4  Means  of  features  for  A01,  E01,  and  F01 
averaged  across  task  loads 


Theta 

Alpha 

Beta 

LL 

1.283 

0.514 

0.866 

ML 

1.237 

0.435 

0.614 

HL 

1.241 

0.487 

0.639 

APPENDIX  B 

FALSE  ALARM  DATA  FROM  CONTROL  CHARTS 
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Table  B.l  False  alarms  of  EWM A-Shewart  control  charts 


VEOG 

Theta 

f2 

Theta 

Fz 

Alpha 

Ft 

Theta 

o2 

Alpha 

Pz 

Alpha 

Average 

A01 

1 

2 

1 

2 

0 

1 

1.167 

E01 

0 

0 

1 

1 

0 

0 

0.333 

F01 

1 

0 

1 

0 

1 

0 

0.500 

Average 

0.667 

0.667 

1.000 

1.000 

0.333 

0.333 

Table  B.2  False  alarms  of  time  series  residuals 
control  charts 

HEOG 

Beta 

F, 

Beta 

Pz 

Beta 

Ts 

Beta 

Average 

A01 

1 

1 

2 

1 

1.250 

E01 

2 

1 

1 

1 

1.250 

F01 

4 

4 

3 

1 

3.000 

Average 

2.333 

2.000 

2.000 

1.000 

Table  B.3  False  alarms  of  MCEMWA 
control  charts 


Fz 

Beta 

Pz 

Beta 

Ts 

Beta 

Average 

A01 

3 

3 

3 

3.000 

E01 

2 

2 

1 

1.667 

F01 

2 

4 

2 

2.667 

Average 

2.333 

3.000 

2.000 
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Table  B.4  False  alarms  of  various  multivariate  control  charts 


Hotelling-T2 

(Select) 

MEWMA 

(Select) 

Hotelling- 

T2  (PCA) 

MEWMA 

(PCA) 

Average 

A01 

4 

2 

4 

3 

3.250 

E01 

4 

1 

2 

2 

2.250 

F01 

6 

4 

4 

2 

4.000 

Average 

4.667 

2.333 

3.333 

2.333 

APPENDIX  C 

PRINCPLE  COMPONENT  FACTOR  PATTERN  ANALYSIS 

Table  C.  1  Principle  component  factor  pattern  of  AO  1 
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pi 

P2 

P3 

P4 

P5 

P6 

P7 

P8 

P9 

P10 

VEOG 

delta 

0.301 

-0.696 

0.295 

-0.036 

0.053 

-0.460 

0.098 

-0.016 

0.017 

-0.026 

theta 

0.353 

-0.634 

0.308 

-0.039 

0.029 

-0.463 

0.236 

0.037 

0.109 

-0.095 

alpha 

0.466 

-0.583 

0.024 

-0.029 

-0.039 

-0.368 

0.359 

0.003 

0.074 

-0.045 

beta 

0.429 

-0.093 

-0.667 

-0.106 

0.136 

-0.049 

0.454 

-0.017 

-0.197 

0.045 

gamma 

0.338 

0.168 

-0.680 

-0.083 

0.118 

-0.002 

0.426 

-0.034 

-0.237 

0.101 

HEOG 

delta 

0.030 

-0.786 

0.129 

-0.129 

-0.074 

-0.082 

-0.314 

-0.024 

0.044 

0.042 

theta 

0.159 

-0.826 

-0.147 

-0.162 

0.019 

0.041 

-0.296 

-0.007 

0.080 

-0.017 

alpha 

0.191 

-0.695 

-0.512 

-0.191 

0.171 

0.202 

-0.161 

0.039 

0.067 

-0.066 

beta 

0.116 

-0.645 

-0.581 

-0.195 

0.271 

0.175 

-0.116 

0.069 

0.075 

-0.081 

gamma 

0.148 

-0.632 

-0.573 

-0.167 

0.290 

0.157 

-0.105 

0.076 

0.071 

-0.071 

N 

U. 

delta 

0.355 

-0.207 

-0.005 

0.497 

0.140 

0.086 

-0.157 

-0.027 

0.036 

0.604 

theta 

0.230 

-0.235 

-0.149 

0.014 

-0.491 

0.021 

-0.288 

-0.324 

-0.568 

-0.028 

alpha 

0.106 

-0.237 

0.250 

-0.191 

-0.120 

0.437 

0.300 

-0.560 

0.235 

0.080 

beta 

0.804 

0.275 

0.027 

-0.103 

0.036 

0.071 

-0.025 

-0.018 

0.063 

0.030 

gamma 

0.808 

0.432 

-0.031 

-0.115 

0.034 

-0.039 

-0.157 

-0.003 

0.058 

-0.006 

r- 

delta 

0.429 

-0.471 

0.114 

0.324 

0.138 

-0.236 

-0.020 

-0.089 

-0.063 

0.360 

theta 

0.398 

-0.637 

0.009 

-0.013 

-0.198 

-0.157 

-0.168 

-0.161 

-0.153 

-0.049 

alpha 

0.450 

-0.356 

0.070 

-0.017 

-0.312 

0.256 

0.198 

-0.305 

0.303 

-0.017 

beta 

0.745 

0.306 

-0.172 

-0.017 

-0.013 

0.095 

0.022 

-0.024 

0.052 

0.094 

gamma 

0.776 

0.337 

-0.096 

-0.040 

0.010 

-0.017 

-0.084 

-0.019 

0.025 

0.135 

Table  C.l — Continued 


delta 

0.378 

-0.155 

0.017 

0.759 

0.167 

0.182 

-0.027 

-0.048 

-0.043 

-0.104 

theta 

0.287 

-0.341 

0.042 

-0.022 

-0.632 

0.140 

0.177 

0.180 

-0.235 

-0.013 

N 

Cm 

alpha 

0.174 

-0.250 

0.590 

-0.313 

0.320 

0.347 

0.115 

0.042 

-0.257 

-0.035 

beta 

0.852 

0.090 

0.213 

-0.110 

0.084 

-0.004 

-0.034 

0.062 

-0.025 

-0.072 

gamma 

0.790 

0.421 

-0.036 

-0.149 

-0.041 

-0.064 

-0.194 

-0.030 

0.070 

-0.065 

delta 

0.338 

-0.120 

0.002 

0.760 

0.183 

0.089 

0.038 

-0.013 

-0.073 

-0.238 

theta 

0.328 

-0.288 

-0.074 

0.223 

-0.470 

0.160 

-0.024 

0.178 

0.171 

-0.181 

H 

alpha 

0.382 

-0. 146 

0.504 

-0.350 

0.234 

0.298 

0.057 

0.020 

-0.152 

-0.019 

beta 

0.845 

0.329 

0.020 

-0.047 

0.011 

-0.079 

-0.010 

0.067 

0.019 

-0.095 

gamma 

0.768 

0.509 

-0.116 

-0.125 

-0.053 

-0.063 

-0.061 

-0.015 

0.003 

-0.068 

delta 

0.301 

-0.087 

0.114 

0.761 

0.222 

0.151 

0.066 

-0.055 

-0.038 

-0.256 

theta 

0.216 

-0.197 

-0.090 

0.195 

-0.517 

0.251 

0.143 

0.504 

0.125 

0.179 

O 

alpha 

0.291 

-0.370 

0.517 

-0.198 

0.214 

0.262 

0.056 

0.335 

-0.155 

0.088 

beta 

0.852 

0.021 

0.213 

-0.040 

0.038 

-0.041 

-0.024 

0.080 

0.007 

0.001 

gamma 

0.827 

0.395 

-0.011 

-0.118 

-0.028 

-0.053 

-0.153 

-0.020 

0.067 

-0.049 
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Table  C.2  Principle  component  factor  pattern  of  E01 


pi 

P2 

P3 

P4 

P5 

P6 

VEOG 

delta 

-0.5189 

-0.6534 

0.144 

-0.2419 

0.0522 

0.2724 

theta 

-0.4243 

-0.7641 

0.0467 

-0.159 

0.1439 

0.2716 

alpha 

-0.433 

-0.6961 

0.1683 

-0.0844 

0.0845 

0.3606 

beta 

-0.6895 

-0.3628 

0.0832 

-0.3572 

0.1927 

0.0261 

gamma 

-0.7851 

0.0523 

0.2799 

-0.2339 

0.272 

-0.0508 

HEOG 

delta 

-0.336 

-0.4305 

0.355 

0.5134 

0.0103 

0.2619 

theta 

-0.3309 

-0.384 

0.4095 

0.6438 

-0.0397 

0.0569 

alpha 

-0.3289 

-0.2771 

0.4425 

0.6549 

0.0255 

-0.2625 

beta 

-0.5119 

-0.1063 

0.4952 

0.3784 

0.2077 

-0.3599 

gamma 

-0.6706 

0.1094 

0.4881 

0.0588 

0.2694 

-0.2786 

N 

ta 

delta 

-0.6035 

-0.4614 

-0.1269 

-0.2895 

-0.1911 

-0.2703 

theta 

-0.3918 

-0.5648 

-0.4679 

-0.114 

-0.076 

-0.3507 

alpha 

-0.4978 

-0.1078 

-0.5945 

0.3151 

0.2686 

-0.0735 

beta 

-0.8873 

0.2272 

-0.1378 

0.0452 

0.1014 

0.0789 

gamma 

-0.9107 

0.2738 

0.1111 

-0.1049 

0.0996 

0.0166 

r- 

delta 

-0.7248 

-0.4215 

0.0103 

-0.2809 

0.0081 

-0.1126 

theta 

-0.4696 

-0.6763 

-0.3183 

-0.1975 

0.0708 

-0.2488 

alpha 

-0.7304 

-0.2067 

-0.4707 

0.1129 

0.1365 

0.0058 

beta 

-0.9062 

0.2418 

-0.0412 

0.0034 

0.0692 

0.104 

gamma 

-0.8862 

0.2518 

0.1427 

-0.0794 

0.1388 

0.0524 
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Table  C.2 — Continued 


delta 

-0.6402 

-0.08 

0.1507 

-0.1932 

-0.5275 

-0.0622 

theta 

-0.4799 

-0.2524 

-0.2622 

0.3264 

-0.4346 

-0.0246 

N 

Oh 

alpha 

-0.4693 

-0.0026 

-0.5046 

0.4817 

0.1778 

0.0789 

beta 

-0.9078 

0.2354 

-0.1689 

0.0688 

0.0452 

0.0739 

gamma 

-0.911 

0.3084 

0.0667 

-0.1007 

0.0427 

0.0329 

delta 

-0.6846 

0.1718 

0.2425 

-0.1408 

-0.3044 

-0.05 1 

theta 

-0.7532 

0.1347 

-0.1896 

0.0684 

-0.1755 

-0.1948 

*Ci 

H 

alpha 

-0.7807 

0.2635 

-0.3097 

0.2122 

0.0648 

0.0046 

beta 

-0.9077 

0.3354 

0.0736 

-0.097 

0.0588 

0.045 

gamma 

-0.8565 

0.3032 

0.189 

-0.1725 

0.0366 

-0.0144 

delta 

-0.6391 

-0.0012 

0.2664 

-0.0898 

-0.5272 

0.0437 

theta 

-0.5422 

0.0021 

-0.1065 

0.3426 

-0.4705 

0.1319 

N 

o 

alpha 

-0.6947 

0.2539 

-0.3355 

0.2498 

-0.0538 

0.1979 

beta 

-0.8691 

0.267 

-0.0338 

-0.0628 

-0.0304 

0.0687 

gamma 

-0.8201 

0.2839 

0.0569 

-0.1676 

-0.008 

0.0435 
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Table  C.3  Principle  component  factor  pattern  of  F01 


PI 

P2 

P3 

P4 

P5 

P6 

VEOG 

delta 

0.6481 

0.4212 

-0.3956 

0.2711 

-0.2277 

0.1072 

theta 

0.6769 

0.4747 

-0.4191 

0.156 

-0.093 

0.166 

alpha 

0.6629 

0.4442 

-0.4431 

0.1988 

-0.1474 

0.1474 

beta 

0.7153 

0.4356 

-0.3604 

0.0943 

-0.1592 

-0.0201 

gamma 

0.7047 

0.3298 

-0.2268 

-0.0884 

-0.1227 

-0.2113 

HEOG 

delta 

0.4053 

0.4785 

-0.007 

-0.2073 

0.5887 

0.0549 

theta 

0.3171 

0.4229 

-0.0119 

-0.2102 

0.7706 

0.1341 

alpha 

0.327 

0.401 

0.0166 

-0.2174 

0.7216 

0.0349 

beta 

0.7429 

0.0939 

-0.0473 

-0.2071 

0.2727 

-0.3319 

gamma 

0.8429 

-0.1172 

-0.0888 

-0.1152 

0.1073 

-0.3839 

N 

U* 

delta 

0.6776 

0.2758 

0.1412 

-0.0529 

-0.2137 

0.1642 

theta 

0.7061 

0.1645 

0.0022 

0.0995 

-0.1612 

0.3521 

alpha 

0.5967 

0.0709 

0.2467 

0.6371 

-0.0346 

-0.0227 

beta 

0.9084 

-0.1756 

-0.129 

-0.0609 

-0.1039 

-0.1901 

gamma 

0.8487 

-0.3597 

-0.1385 

-0.1057 

-0.0291 

-0.2219 

r- 

delta 

0.792 

0.0362 

0.0809 

-0.0665 

-0.0444 

-0.198 

theta 

0.8477 

0.1565 

-0.18 

0.1385 

-0.0355 

0.0885 

alpha 

0.7383 

0.1694 

-0.0532 

0.4833 

0.0014 

-0.0765 

beta 

0.8639 

-0.11 

-0.2735 

-0.0184 

-0.0179 

-0.0425 

gamma 

0.8661 

-0.2598 

-0.2619 

-0.0743 

-0.0073 

-0.078 
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Table  C.3 — Continued 


delta 

0.5817 

0.3189 

0.4502 

-0.3692 

-0.3161 

-0.0487 

theta 

0.5812 

-0.005 

0.3404 

-0.0409 

0.0119 

0.2805 

N 

alpha 

0.3621 

-0.1135 

0.4305 

0.6667 

0.2255 

-0.1346 

beta 

0.8372 

-0.4329 

-0.0042 

-0.0682 

0.1052 

-0.1091 

gamma 

0.7565 

-0.5424 

-0.0852 

-0.0934 

0.0703 

-0.0969 

delta 

0.5121 

0.3098 

0.5481 

-0.4078 

-0.2885 

-0.1013 

theta 

0.6887 

-0.0558 

0.4622 

-0.1293 

-0.0936 

0.0457 

«/•> 

H 

alpha 

0.4858 

-0.068 

0.5637 

0.483 

0.1444 

-0.1009 

beta 

0.6636 

-0.559 

-0.0096 

-0.1271 

0.1269 

0.1912 

gamma 

0.5974 

-0.6202 

-0.1568 

-0.1055 

0.1261 

0.2266 

delta 

0.4669 

0.3024 

0.5679 

-0.383 

-0.2818 

-0.1089 

theta 

0.3816 

0.222 

0.4075 

-0.2627 

-0.013 

0.3271 

fS 

C 

alpha 

0.2958 

-0.0041 

0.5662 

0.5619 

0.1806 

0.0296 

beta 

0.7242 

-0.4314 

0.0924 

-0.0718 

0.0529 

0.3149 

gamma 

0.6222 

-0.5409 

-0.0505 

-0.0899 

-0.0604 

0.3376 
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